
## This file provides some commonly used helper functions


## standardize inputs
## gelman=TRUE divides by 2 standard deviations

zstd <- function(x, gelman=FALSE) {
	if (gelman==FALSE){
		if (is.vector(x)==TRUE) {
			out <- (x - mean(x, na.rm=TRUE))/sd(x, na.rm=TRUE)
		}
		else {
			zstd.fun <- function(x) {(x - mean(x, na.rm=TRUE))/sd(x, na.rm=TRUE)}
			out <- apply(x, 2, zstd.fun)
		}
		return(out)
	}
	if (gelman==TRUE){
		if (is.vector(x)==TRUE) {
			out <- (x - mean(x, na.rm=TRUE))/(2*sd(x, na.rm=TRUE))
		}
		else {
			zstd.fun <- function(x) {(x - mean(x, na.rm=TRUE))/(2*sd(x, na.rm=TRUE))}
			out <- apply(x, 2, zstd.fun)
		}
		return(out)
	}

}


## center variables
## group=TRUE does group specific centering

center <- function(x, group=FALSE, groupid=NULL) {
	if (group==FALSE){
		if (is.vector(x)==TRUE) {
			out <- x - mean(x, na.rm=TRUE)
		}
		else {
			center.fun <- function(x) {x - mean(x, na.rm=TRUE)}
			out <- apply(x, 2, center.fun)
		}
	}
	else if (group==TRUE){
		if (is.vector(x)==TRUE) {
			out <- rep(NA, length(x))
			for (i in unique(groupid)){
				out[groupid==i] <- x[groupid==i] - mean(x[groupid==i], na.rm=TRUE)
			}
		}
		else {
			out <- matrix(NA, ncol=ncol(x), nrow=nrow(x))
			center.fun <- function(x) {x - mean(x, na.rm=TRUE)}
			for (i in unique(groupid)){
				out[groupid==i,] <- apply(x[groupid==i,], 2, center.fun)
			}
		}
	}
	return(out)
}


## function to generate data list for JAGS
gen.list <- function (vector.of.variable.names) {  
	sapply(vector.of.variable.names, function(x) get(x))
}




#  sum.coda.R
#  calculates posterior summaries of coda samples
#
#  D.S., Mannheim 2011-08-14.
# 26/12/11: print HPDs by default, say CI=TRUE for cridible int.
# 

		
sum.coda<- function (obj, start=1, end=NA, conf = 0.95, ci=FALSE, digits = 3) 
{
	hpd <- function (x, hpdlvl=conf) {
    n <- length(x)
    m <- max(1, ceiling((1-hpdlvl) * n))
    y <- sort(x)
    a <- y[1:m]
    b <- y[(n - m + 1):n]
    i <- order(b - a)[1]
    structure(c(a[i], b[i]))
	}
    
	if(is.na(end)==TRUE){
		end <- dim(obj)[1]
		cat("Chain length:", end, "(Default start and end used)", "\n\n")
		}
	else if(is.na(end)==FALSE){
		cat("Chain length:", end-start, "\n\n")
		}
	
 	 means <- apply(as.matrix(obj[start:end, ]), 2, mean)
 	sds <- apply(as.matrix(obj[start:end, ]), 2, sd)
        pdens <- apply(as.matrix(obj[start:end, ]), 2, hpd)
        pdens2 <- apply(as.matrix(obj[start:end, ]), 2, quantile, prob=c((1-conf)/2, 1-(1-conf)/2))
  
    if (ci==TRUE){
		mat <- round(cbind(means, sds, pdens2[1, ], pdens2[2, ]), digits)
		colnames(mat) <- c("Mean", "    sd", "  low", "high")
		}
	if (ci==FALSE){
		mat <- round(cbind(means, sds, pdens[1, ], pdens[2, ]), digits)
		colnames(mat) <- c("   Mean", "    sd", "  low", "high")		
	}
    return(mat)
}


